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Abstract. The domain integral equation method with its FFT-based matrix- vector products is 
a viable alternative to local methods in free-space scattering problems. However, it often suffers from 
the extremely slow convergence of iterative methods, especially in the transverse electric (TE) case 
with large or negative permittivity. We identify the nontrivial essential spectrum of the pertaining 
integral operator as partly responsible for this behavior, and the main reason why a normally efficient 
deflating preconditioner does not work. We solve this problem by applying an explicit multiplicative 
regularizing operator, which transforms the system to the form 'identity plus compact', yet allows 
the resulting matrix-vector products to be carried out at the FFT speed. Such a regularized system 
is then further preconditioned by deflating an apparently stable set of eigenvalues with largest mag- 
nitudes, which results in a robust acceleration of the restarted GMRES under constraint memory 
conditions. 
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1. Introduction. There is a steady interest in the numerical simulation of the 
electromagnetic field in inhomogeneous media. The methods can be roughly divided 
into two categories - local and global - in accordance with the governing equations. 
The local methods based on the differential Maxwell's equations [9j [101 HH [H] are 
generally more popular due to the sparse nature of their matrices and the ease of pro- 
gramming. Alongside the local methods, global methods based on an equivalent inte- 
gral equation formulation are also often employed, especially in free-space scattering. 
If an object is large and homogeneous or has a perfectly conducting boundary then the 
problem is usually reduced to a boundary integral equation with the fields (currents) 
at the interfaces being the fundamental unknowns |13) . If an object is continuously 
inhomogeneous or is a composite consisting of many small different parts, the most 
appropriate global method is the domain integral equation (DIE), in two dimensions, 
or the volume integral equation, in three dimensions [TB I H5 1 IT9 ] I2T1 \22\ [25 ] 127 ] [28 ] . 
Although the global methods produce dense matrices, they are generally more stable 
with respect to discretization than the local ones, and the convolution-type integral 
operators sometimes allow to compute matrix-vector products at the FFT speed. 
These properties make the DIE method a viable alternative to local methods for 
certain free-space scattering problems. 

The main difficulties with the DIE method are the non-normality of both the 
operator and the resulting system matrix, inherent to frequency-domain electromag- 
netic scattering, and the extremely slow convergence of the few iterative methods that 
can be applied with such matrices. Numerical experiments consistently show that the 
GMRES algorithm is the best choice [11] , and it can be proven that the full (un- 
restarted) GMRES will evetually converge with any physically meaningfull scatterer 
and incident field [23] . Most of the practically interesting problems, however, involve 



'School of Electrical and Computer Engineering, National Technical University of Athens, Athens 
15773, Greece (zourosamail.ntua.gr). 

^Numerical Analysis, DIAM, Faculty of Electrical Engineering, Mathematics, and Computer Sci- 
ence, Delft University of Technology, 2628 CD Delft, The Netherlands (n. v. budkoStudelft.nl). 



2 



G. P. ZOUROS AND N. V. BUDKO 



the number of unknowns prohibitive for the full GMRES, and its restarted version 
is used. Unfortunately, sometimes, especially in TE scattering with large permittiv- 
ities/dimensions, restarted GMRES as well Bi-CGSTAB converge very slowly or do 
not converge at all. 

Several successful preconditioning strategies have been proposed for local meth- 
ods, see e.g. [SJ [III]) and boundary integral formulations [TJ [S] [SJ [Tj [55]. Whereas 
the domain integral equation method is lagging behind in this respect [SJ [TTJ [T2] HH] ; 
reporting convergence and acceleration thereof mainly for low-contrast and/or small 
objects. One of the difficulties in designing a suitable multiplicative preconditioner 
for this method is that it needs to be cither sparse or have a block- Toeplitz form 
to be able to compete with local methods in terms of memory and speed. Recently 
a deflation-based preconditioner has been proposed in [24] for accelerating the DIE 
solver in the TM case. A thorough spectral analysis of the system matrix reported 
in [21] showed that the outlying eigenvalues (largest in magnitude) are responsible 
for a period of stagnation of the full GMRES, and that their deflation accelerates the 
iterative process. As we show here this strategy, however, does not work in the TE 
case. Moreover, using the deflation on top of a full GMRES puts an even higher strain 
on the memory. 

To further extend the applicability of the DIE method towards objects with higher 
contrasts/sizes, in this paper we propose a two-stage preconditioner based on the 
regularization of the pertaining singular integral operator and subsequent deflation of 
the largest eigenvalues. To derive a regularizer we employ the symbol calculus, which 
also provides us with the essential spectrum of the original operator. Our regularizer 
may be viewed as a generalization of the Calderon identity recently employed for 
preconditioning of the boundary integral equation method [TJ [2[ [6] . As an illustration 
and a proof of the concept we compute the complete spectrum of the system matrix for 
a few typical objects of resonant size. We analyze the difference in spectra between 
the TE and TM polarizations and demonstrate that the discretized version of the 
regularizer indeed contracts a very dense group of eigenvalues making the regularized 
system amenable to deflation. Unlike the previously mentioned boundary integral 
equation approach pQ the DIE method permits a straightforward discretization of 
the continuous regularizer without any additional precautions. As our regularizer is 
a convolution-type integral operator, all matrix-vector products can be carried out 
with the FFT algorithm. The subsequent deflation of the largest eigenvalues can be 
achieved with the standard Matlab's eigs routine. Most importantly, these largest 
eigenvalues and the assiociated eigenvectors turn out to be rather stable, so that one 
can considerably accelerate the eigs algorithm by setting a modest tolerance. We have 
performed a number of numerical experiments which have consistently demonstrated 
that given a fixed amount of computer memory it is better to spend some part of it on 
deflation than all of it on maximizing the dimension of the inner Krylov subspace of 
the restarted GMRES. Keeping the restart parameter of the GMRES in the order 
or larger than the dimension of the deflation subspace gives a speed-up which grows 
with the object size and permittivity. 

Although the DIE method (both TM and TE versions of it) has been around 
in engineering community for many years [TU [13 HH [19l EU US [25l [27l [28], the 
full theoretical analysis of the two-dimensional case (in the strongly singular form) 
is still missing. The reason might be the presence of special functions in the kernel 
of the integral operator, which makes calculations more tedious than in the three- 
dimensional case [UdS]. Therefore, before going into the numerical details, we devote 
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some space here to the analysis of the two-dimensional singular integral operator. As 
a result we show an almost complete equivalence of the TE scattering to the more 
general three-dimensional scattering problem. Explicit numerical calculations of the 
matrix spectrum for resonant scatterers presented here are generally not possible with 
three-dimensional objects of that size. Whereas, the established spectral equivalence 
indicates that the proposed preconditioner will work with three-dimensional problems 
as well. 

2. Domain integral equation. We consider the two-dimensional frequency- 
domain electromagnetic scattering problem on an inhomogeneous object of finite spa- 
tial support exhibiting contrast in both electric and magnetic properties. We assume 
that the object (a cylinder) as well the source of the incident field (a plane wave, 
current carrying line, etc.) are invariant in the a^-direction of the spatial Carte- 
sian system. The Maxwell equations describing the total electromagnetic field in the 
(xi, X2)-plane in the presence of an inhomogeneous magneto-dielectric scatterer can 
be written as 



-d 2 E 1 (x,u>) - 
d 2 H 1 (x,ui) 



-d 2 H 3 (x,uj) 
diH 3 (x,u>) 
d 1 E 2 (x,uj) - 

- diH 2 {x,u) 



- iu>e(x,oj)Ei(x,ui 

- iuie(x, u>)E 2 (x, lj 

iui/J,(x, Ul)H 3 (x, UI 

- iuie{x, ui)E 3 (x, us 
iw/x(x, uj)Hi(x, ui 



—d\E 3 {x,u>) — iu)n(x,uj)H 2 (x,u) 



= -Ji(x,w), 
= -J 2 (x ) ui) ) 
= 0, 

= -J 3 (x,w), 
= 0, 
= 0, 



(2.1) 



(2.2) 



where we introduce the two-dimensional position vector x = (xi, x 2 ); J is the electric 
current density; e and fj, are the possibly complex-valued dielectric permittivity and 
magnetic permeability; and ui is the angular frequency. It is easy to notice that the 
first three equations, (|2.ip . are decoupled from the last three, (|2.2[) . Reflecting the 
fact that the electric field strength is transverse to the axis of invariance, equations 
(12. 1|1 are said to describe the transverse electric or TE case (also somewhat confusingly 
called H-wave polarization) . Analogously, the three equations (|2.2j) correspond to the 
transverse magnetic, or TM case (also known as E-wave polarization). 

The integral formulation of the scattering problem is obtained by first introducing 
the so-called incident field (E ln , H ln ), which has the same physical source as the total 
field, but satisfies some simpler version of the above equations. Usually, the incident 
field is considered in a homogeneous isotropic background medium with parameters 
£b(w) and /ib( w )- Then, it is fairly easy to see that the scattered field, defined as the 
difference between the total and the incident fields 



E sc (x,w) = E(x,w) - E in (x,w), H sc (x,u;) = H(x,w) - H ln (x,w), 
also satisfies Maxwell's equations in the same homogeneous medium, i.e., 





-d2H s 3 c (x) 


- iu)S h El c (x) 


= - ^i nd (x), 




diH s 3 c (x) 


- ioje h El°(x) 


- -4 nd (*), 


d 2 ET(x)- 




- iuifibH 3 c (x) 


= -*T d (x), 


d 2 Hl c (x) 


-dxH s 2 c (x) 


— iuiebE 3 ic (x) 


= -4 nd W, 




d 2 E?(x) - 


- iuj[j, h Hl c (x) 


= -Al nd (x), 




-9xE^(x) 


- iujfibH 2 c (x) 


= -*T d (x), 



(2.3) 



(2.4) 



(2.5) 
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where the implicit dependence on ui has been omitted for simplicity. The physical 
sources of the scattered field are the induced current densities, also known as contrast 
currents: 



4 nd (x) = -tw [e(x) - e b ] £? fc (x), k = 1, 2, 3; 
iC d (x) = -iu [/x(x) - Mb ] H m (x), m = 1,2, 3. 

It is convenient to write the Maxwell equations in the following matrix form: 



(2.6) 



— iuJEh 





-82 ' 




'Ef' 




~-Jf d ~ 





— iuJSh 


di 




Ef 




rind 


_ -d 2 


d x 


—iuj/Jh 




Hf_ 




-Kf\ 


— iw/Xb 





d 2 




~Hf~ 











-d x 




Hf 




-Kf d 


d 2 


-d x 


—iijje-x, 




_Ef_ 




rind 
_~ J 3 _ 



(2.7) 



(2.8) 



The solution of these equations can be obtained by first transforming l|2.7[) - (j2.8[) to 
the k-domain with the help of the two-dimensional spatial Fourier transform, then 
solving the resulting algebraic systems, and finally transforming the result back to 
the x-domain. The obtained solution expresses the scattered field in terms of the 
induced currents and, via (|2.6j) . in terms of the unknown total field. Subsequently, 
using (|2.3[) . one can replace the scattered field with the difference between the total 
and the incident fields, thus arriving at the following integro-differential equations 
with the total field as a fundamental unknown: 



(2.9) 



'Ef' 




'Ef 




kl + dl dxd 2 


-iLo^ h (-d 2 ) 




9* {XeEl) 


Ef 




Ei 




d 2 dx kl + di 


—iu!/j,hdx 




9* (XcE 2 ) 


Hf 




H 3 




-iuJEb(-d 2 ) -iujebdi 


k 2 




9 * (XmH 3 )_ 


'Hf 




~H{ 




' kl + dl dxd 2 


-iujehd 2 




9 * (XmHx) 


Hf 




H 2 




d 2 dx kl + di 


-iuis h {-dx) 
k 2 

K b 




9 * (XmH 2 ) 


Ef 




E 3 




-iaj[ibd 2 —iu)^b(—dx) 




_9*(XeE 3 )_ 



(2.10) 



Here we have introduced the wavenumber of the homogeneous background medium 
kl = Lu 2 EbHb and the normalized electric and magnetic contrast functions: 



Xc(x,w) 
Xm(x,w) 



e(x,a;) 

£b 
/x(x, Ul) 

Mb 



1. 



1. 



(2.11) 
(2.12) 



The star (*) denotes the two-dimensional convolution, e.g., 



g * ( Xe Ex) = f ff (x - x')Xe(x')^i(x')dx', (2.13) 
The scalar Green's function of the two-dimensional Helmholtz equation is given by 

g(x 1 u>) = ±H< i 1 \k b \x.\), (2.14) 



where Hq is the zero's order Hankel function of the first kind. 
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3. Operator, symbol, spectrum, and regularizer. Although the integro- 
differential equations (|2.9p - (|2.10l) can to a certain extent be analyzed directly on the 
pertaining Sobolev spaces, we prefer to transform them into singular integral equations 
and work on the Hilbert space L\ (D) of vector- valued functions with spatial support 
on D. 



As the two problems (|2.9|) and (|2.10l) are mathematically identical, up to some 
constants and non-essential sign changes in the differential matrix operator, we shall 
concentrate on just one of them, say the TE case (|2.9[) . A standard singular integral 
equation is obtained from (|2.9I) by carrying out the spatial derivatives of the weakly 
singular integrals (|2.13[) , and writing the result as 



Du + GXu + KXu, 



(3.1) 



where D and X are "diagonal" multiplication operators, 67 is a principal- value singular 
integral operator, and K is a compact integral operator. In our case, things become a 
little more complicated due to the matrix- valued kernel and the involvement of special 
functions, but eventually one can put equation (12.91) in the standard form as 



Ef 

Hip 



S 
5 
0/ 



E 2 
H 3 



p. v. 



Gu Gx2 
G21 G 2 2 




Kxx K12 Ki 3 
K 2 i K 2 2 K23 
K31 K32 K 33 



X e Ei 
X C E 2 
X m H 3 

X e Ex 
X C E 2 
X m H 3 



(3.2) 



Here, the first matrix operator on the right contains the identity operator /, and two 
operators S, which denote pointwise multiplication with the following function: 



s ( x ) = 1 + 



(3.3) 



Operators X e / m denote pointwise multiplication with the corresponding contrast func- 
tions Xc/m( x )- The second term on the right is the principal value of a convolution- 
type integral operator (with circular exclusion area around the singularity, which 
tends to zero in the limit), where the nonzero elements of the matrix- valued kernel 
are 



Grim ( x ) 



271-1x1 



■ [29 n 9 m - S nm ] , n,m= 1, 2; 



(3.4) 



with 9„ 



|x| being a Cartesian component of a unit vector. The last term on the 
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right in (13. 2[) is also a convolution with the matrix elements of the kernel given by 

1 



2tt|x| 2 
4 



b ^ 1) (A^|x|)[M n 



<5nm] , n,m = 1, 2; 



#13 = - 



31 



A" 



23 



^^( fcb |x| 
W£ bfcb% if (i) (fcb|x|) 

^(fcblxl), 



4 

4 

web^b^ 1 



33 



-^i (1) (fcb|x| 



^^(fcblxl). 



ife 2 
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(3.5) 

(3.6) 
(3.7) 
(3.8) 
(3.9) 
(3.10) 



As can be deduced from asymptotic expansions of special functions, all K nm 
kernels are weakly singular, i.e., their singularity is less than that of the factor |x — 
x'|~ 2 , and the corresponding integrals exist in the usual sense. It is also clear that the 
matrix integral operator KX is compact on L 2 (D) . The remaining strong singularity 
is explicitly contained within the G nm kernel. The latter kernel also features a very 
special tensor 29 n 9 m — 8 nm , which guarantees the existence of the integral in the sense 
of principal value. The standard assumption, sufficient to guarantee the boundedness 
of the matrix singular integral operator G and the equivalence of the integral equation 
(|3.ip to the original Maxwell's equations, is that of Holder continuity of the dielectric 
permittivity and of all components of the incident field vector for all x £ M 2 . Although 
in practice the equation is posed and solved on L\ (D) only, the existence theory was 
developed for domains without an edge [20], i.e., on (R 2 ). Samokhin, however, 
has shown, [23] , that as far as the question of existence is concerned all the operators 
can be extended to L^(R 2 ) without loosing the important compactness of weakly 
singular integral operators. 

At the core of the Mikhlin-Prossdorf approach, sucessfully applied in [U [23] to the 
three-dimensional volume integral equation, is the symbol calculus. It allows to study 
the existence of a solution and shows a way to reduce the original singular integral 
operator to a manifestly Fredholm operator of the form 'identity plus compact'. The 
symbol of an integral operator in the present two-dimensional case is a matrix-valued 
function $ nm (x, k), n,m = 1,2,3, x G R 2 , and k e M 2 . Construction of the symbol 
follows a few simple rules explained in [301 1231 H] : the symbol of the sum/product of 
operators is the sum/product of symbols; the symbol of a compact operator is zero; 
the symbol of a multiplication operator is the multiplier function itself; the symbol of 
the strongly singular integral operator is the principal value of the Fourier transform 
of its kernel (that is where the second argument k - the dual of x - comes from). 
It follows, in particular, that the symbol of the matrix operator is a matrix-valued 
function. Applying these rules to (13.21) we obtain 



*(x,k) = 



S (x) 






S (x) 




Gu(k) 
G 2 i(k) 




Gi 2 (k) 
G 22 (k) 




Xe(x) 








Xe(x) 







Xm(x) 



(3-11) 
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The Fourier transforms G nm (k) of the kernels G Tl 
dimensional 'spherical' functions of order p: 



t (x) are found as a series of two- 



G nm (k) 



p=0 



a (2)y(2) 
u p 1 p,2 



(3.12) 



which, in general, results in a different expression for each combination of n and m. 
The particular form of (|3. 12|) reflects the fact that G nm is a function of the direction 
of k only and it does not depend on its magnitude (see [U [20] for a formal proof of this 
fact). In two dimensions, with (f> denoting the directional angle of the unit vector k/|k| 
we have: Y$ = sin(p0), = cos(p0), and 72 , p = m p T(p/2)/T((p + 2)/2). The 
expansion coefficients in (|3.12l) can be easily recovered by representing the original 
kernels G„ m (x) in the following form: 



Gnm (x^) — . |2 fnm 



If we now express the characteristic f nm as 

^ nm ( Ixl ) = ^ sin ^ + °P 2) cos (^) 



(3.13) 



(3.14) 



where 4> is the directional angle of the unit vector x/|x|, then the expansion coefficients 
in (|3.12[) are those of A3. 14)) . Again we note that for our matrix- valued kernel this has 
to be done for each combination of n and m separately. Eventually we arrive at 



Gn(k) 
G 2 i(k) 




G 12 (k) 
G 22 (k) 




cos 2 (</>) - 1/2 
sin(0) cos(4>) 




sin(0) cos(</>) 
sin 2 (<?<>)- 1/2 




(3.15) 



"1 





0" 




cos 2 (0) sin(< 


t>) cos(^) 


0" 







1 





+ Xc(x) 


sin(^) cos(0) si 


n 2 (0) 





= I + Xo(x 








1 
















Finally the symbol of the complete operator is 



*(x,< 



(3.16) 

where matrix Q is a projector, i.e., Q 2 = <Q>. 

Several important conclusions can be made already at this stage. First of all 
we notice an almost complete equivalence of the symbol in the two-dimensional TE 
case with the previously obtained symbol of the three-dimensional scattering operator 
[H [53] . The necessary and sufficient condition for the existence of a solution is now 
readily obtained as the condition on the invertibility of 4>, which is 



e(x, oj) 7^ 0, x e 



(3.17) 



In the TM case (|2.10[) the symbol can be derived along the same lines and turns 
out to be I + Xm(x)Q, with the existence condition being /x(x, w) 7^ 0, x G R 2 . It 
is interesting to note the difference in these existence conditions with respect to the 
three-dimensional case, where both e(x) and /Lt(x) are required to be nonzero at the 
same time. 
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The explicit inverse of the symbol matrix is found to be 

*- 1 (x^)=I + x ;(x)Q, (3.18) 

where Xe( x ; w ) — £ b/e(x,w) — 1. The inverse of the symbol is the symbol of the 
regularizer - an operator, which reduces the original singular integral operator to the 
form 'identity plus compact'. Since has the same form as we conclude that our 
original operator with its electric contrast function changed into x' e m ay be employed 
as a regularizer. Thus, if A(x c ) is the original operator from (|3. ip — (|3.2|) . and A(x' e ) 
denotes the same operator with the contrast function replaced by x'm then we have 

A( X ' e )A(Xe) = I + K, (3.19) 

where if is a generic compact operator. This looks strikingly similar to the Calderon 
identity [U [H [B] , and can be viewed as an extension of the latter to the case of an 
inhomogeneous penetrable scatterer. 

Regularizer eliminates the essential spectrum of A(x c )- Indeed, even if A(x c ) has 
a nontrivial essential spectrum, densely spread over a part of the complex plane, the 
spectrum of I + K consists of isolated points only. As can be deduced from the symbol 
(|3.16p , the essential spectrum of the operator A(xc) m the TE case is given by 

A css = ^, xel 2 , (3.20) 
£b 

and by A oss = /j,(x)//j,b, x € I 2 in the TM case. From the condition of the uniqueness 
of the solution we also obtain a traditional [18] wedge-shaped bound on the possible 
discrete eigenvalues: 



Ime(x) — [Im e(x) + Im £b] Re A 
+ [Ree(x) -ReebJImA + ImeblA^ < 0, xefl 



(3-21) 



4. GMRES convergence and matrix eigenvalues. Let us introduce a uni- 
form two-dimensional TV- node grid over a rectangular computational domain D. For 
simplicity let the grid step h be the same in both directions. Typically the grid step is 
chosen in accordance with the following empirical rule: find out the highest local value 
of the refractive index n = max{ Y / e'(x)/x'(x)]/eb/ib}, x G D, where prime denotes 
the real part; choose an integer k - the number of points per smallest wavelength; 
and, finally, compute the grid step as h = Ab/ (kn), where Ab is the wavelength in the 
background medium. We set our discretization here at the traditional level of k = 15 
points, which, as we verified, gives a small global error (in the order of 10 -5 ) with 
respect to the analytical solution for a homogeneous circular cylinder. For rectangu- 
lar grid-conforming boundaries the error with this discretization rule could be even 
smaller, since most of the error for a circular cylinder comes from a poor geometrical 
representation of the boundary and/or bad approximation of the area of the circular 
cross-section. 

Applying a simple collocation technique with the mid-point rule we end up with 
an algebraic system of the order 3N x 37V: 

Au = 6, (4.1) 



with the following structure: 



All 


An 


A13 




Ml 






A 2 i 


A 22 


^23 








62 


A 3 i 


^32 


^33 




U3 




63 



(4.2) 
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where the matrix elements are 

[Au] nm = - &b^ 2 Xc(x m ) 
+ 



ri in 

i 



fit 



nm^t.nm 



Hi ^ \k\yTnm) 



(4.3) 



1, 2, m =/= n 



\A, 



iirk h h (1) r 



Xc( x n)j 



1,2 



(4.4) 



„, -f^l ^(fcb?"nm) . Hq '(fcb^nm) 



(4.5) 



[^4l3]nm = !W/Jb/i 2 Xm (x m ) ^fff (^b^ml^^m (1 ~ $ n m) 



[A 3 2] nm = -iu£bh 2 Xc(x m ) l -^H[ 1 \k h r nm )9 1 ^ nm (l - S nrn ) 



(4.6) 

(4.7) 
(4.8) 



[A 



23 inn — ~ Kl32 \nm 



[A, 



(4.9) 



(4.10) 



[^4-33]rm — 1 



1 - 



ink^h 



H^(k h h/^) 



Xm( x n) 



In these formulas: r r , 



(4.11) 



1, 2, n, m 



1,...,N. Here X£ yU denotes the Cartesian component of the two-dimensional po- 
sition vector x ra pointing at the nth node of the grid. The unkhown vector u = 
[ui, u 2 , u 3 ] T contains the grid values of the TE total field [i?i(x m ), _E 2 (x m ), iJ 3 (x TO )] T , 
m = 1, . . . , N. The right-hand side vector b = [bi, b 2 , &3] 71 contains the grid values of 
the incident field. For example, the TE plane wave impinging to the object at an angle 
tp with respect to xi-axis is represented by: [bi] n = exp[ifcb(£i,n cos ip + %2,n sin^>)], 
[&a]n = -[6i]n^bSin^/(w£b), and [6 3 ]„ = [&i]„fc b cos ip/(we h ), n=l,...,N. 

Both the system matrix and the original integral operator are neither symmetric 
nor normal. Therefore, the range of iterative methods applicable to the problem is 
extremely limited, with (restarted) GMRES typically showing the best performance. 
Figure |4~T1 shows the GMRES convergence histories on the original system Au = b for 
three physically distinct scatterers (all with Xm = 0) and two different polarizations, 
corresponding to the TE and TM cases. The three scatterers are: a lossless object 
with large positive real permittivity, e/eh = 16, an object with small losses and a 
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Iterations Iterations 

Fig. 4.1. Convergence of restarted GMRES applied to the original system with (dashed) and 
without (solid) deflation of the largest eigenvalues for three physically distinct scatterers. Deflation 
accelerates GMRES in the TM case (right), but does not work in the TE case (left), lossless 
object - red lines, object with small losses and negative real part of permittivity - green lines, and 
inhomogeneous object - blue lines. 



large negative real part of the permittivity, e/eb = — 16 + zl.5, and an inhomogeneous 
object consisting of two concentric layers with the inner part having large losses, 
e/sb = 2.5 + i20, and the outer layer having a large positive real permittivity, e/sb = 
16. All objects are cylinders with square cross-sections with the side length a = At,. 
The inner core of the inhomogeneous object has also a square cross-section with the 
side length a/2. The restart parameter of the GMRES algorithm is equal to 40 
(we shall study the influence of this parameter in more detail later). In Figure POl 
(left) we illuminate the scatterers with a plane wave whose electric field vector is 
orthogonal to the axis of the cylinder (TE) , whereas in Figure I4.1f right) the electric 
field vector of the incident plane wave is parallel to the axis of the cylinder (TM). Since 
all three objects do not exhibit any magnetic contrast with respect to the background 
medium, the integral operators are different for these TE and TM cases - see (|2.9j) . 
(1 2 . 1 [) . The latter contains now only the weakly singular part and is, therefore, of the 
form 'identity plus compact'. 

Next to the convergence plots of the original system in Figure 14.11 we present the 
convergence histories of the restarted GMRES with a deflation-based right precondi- 
tioner (dashed lines). As was suggested in [21], deflating the largest eigenvalues of A 
in the TM case may significantly accelerate the convergence of the full (un-restarted) 
GMRES. The preconditioner deflating r largest eigenvalues of A is constructed by 
first running the eigs algorithm to retrieve these r eigenvalues and the corresponding 
eigenvectors. With the help of the modified Gramm-Schmidt algorithm we further 
build an orthonormal basis for the retrieved set of eigenvectors. Let this basis be 
stored in a 3iV x r matrix V. Then, the right preconditioner [8] has the form 

P- 1 =h N + V[T- 1 -I r ] V*, (4.12) 

with 

T = V*AV. (4.13) 

As was observed in |24j for the TM case, although the system matrix is not normal, 
the matrix V has full column rank, and the matrix T has a decent condition number. 
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Table 4.1 

Restarted GMRES with fair memory usage for three different test scatterers, TM case. Original 
and deflated systems. 



System 


e/e h 


r, deflated 
eigenvalues 


restart 


CPU time, 
seconds 


Iterations, 
tol. lO" 8 


Speed- 
up 


A 


16 





40 


168 


7311 


1.0 


AP- 1 


16 


28 


12 


78 


4230 


2.2 


A 


-16 + il.5 





40 


0.7 


27 


1.0 


AP- V 


-16 + il.5 


28 


12 


0.5 


14 


1.4 


A 


16 and 2.5 + i20 





40 


23 


1000 


1.0 


AP- 1 


16 and 2.5 + i20 


28 


12 


14 


749 


1.6 



15 
5 

-5 
-15 
-25 



-35 



-25 



15 
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-5 
-15 
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(a) TE case. 





















t 


\ ... 





-25 -15 -5 5 15 25 

(b) TM case. 
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Fig. 4.2. Eigenvalues of the system matrix A for three physically distinct scatterers in the 
TE (top) and TM (bottom) cases. Horizontal axis - real part, vertical axis - imaginary part. Also 
shown ( squares, triangles and circles ) the values of the relative permittivity, including the background 
medium, in each scattering configuration. 



We confirm this behavior in the present TE case. The term 'deflation' refers to 
the fact that the spectrum of AP' 1 looks like the spectrum of A with r largest 
eigenvalues moved to the point 1 + i0 of the complex plane, and the corresponding 
eigen-subspace has been projected out. Here, even with the restarted GMRES, we 
see some improvement in the TM case Figure ICT right) . The TE case, however, is 
not affected by deflation (it may even become worse). 

Table 14.11 and Table 15.11 summarize the numerical results corresponding to the 
test objects in the TM and TE cases. The 3 AT x r matrix V needs to be stored in the 
memory. For the sake of fair comparison, we divide our memory between the inner 
Krylov subspace of the restarted GMRES and V, i.e., if we set restart = k, k > r 
with the original matrix A, then we use restart = k — r with the preconditioned 
matrix AP' 1 . 
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Fig. 5.1. Eigenvalues of the regularized matrix A-^A, corresponding to the spectra of A in the 
TE case, see Figure [4-^\ (top). Red circles show the estimated outer bounds of the dense cluster - 
deflation bounds. 



One can understand the peculiar lack of progress in the TE case by consider- 
ing Figure 14.21 which shows the spectra of the system matrix for the TE and TM 
cases corresponding to the numerical experiments of Figure 14.11 As can be seen, the 
difference between the spectra is in the presence of dence line segments [l,e(x)/eb], 
x € R 2 and some extra isolated eigenvalues inside a distribution remotely resembling a 
shifted circumference (the latter is much more pronounced with homogeneous objects 
at higher frequencies, see [23] )• The dense line segments in the TE spectrum contain 
a very large number of closely spaced eigenvalues, and the deflation process simply 
stagnates when it reaches the outer ends of these lines, corresponding to e(x)/eb, 
x e D. Indeed, in the TE case we have seen no change in the convergence with AP~ X 
when we tried to deflate any eigenvalues beyond those that were present outside a 
circle with the center at zero and the radius reaching the largest value of |e(x)/eh|, 
x £ R 2 for each particular scattering configuration. At the same time in the TM case 
we observed consistent improvement in convergence when deflating all eigenvalues 
outside the unit circle, see Table |4~T1 

5. Preconditioning. In the three-dimensional case it was conjectured [3] that 
the dense line segments, which we now observe in the TE spectrum as well, are the 
discrete image of the essential spectrum, A css = e(x)/eb, of the corresponding singular 
integral operator. The absence of these lines in the TM spectrum is just another 
indirect confirmation. Yet, we do not have a formal proof, since there is no such thing 
as essential spectrum with matrices, and its role and transformation in the process of 
discretization is unclear. Anyway, the process of deflation is obviously hindered by 
those dense line segments and, assuming that they are due to the essential spectrum, 
we already know how to 'deflate' them. On the continuous level this can be achieved 
by regularization. We shall, therefore, apply the discretized regularizer ^4r, i.e., a 
discrete version of the operator A(x' c ), so that the problem becomes ArAu = Ar&. 
Notice that no additional storage is required, since the original Green's tensor can 
be re-used. Matrix-vector products, however, do become twice more expensive to 
compute. 

To see how this procedure affects the matrix spectrum we compute the eigenvalues 
of ArA corresponding to the matrices analyzed in Figure l4~2l (top). The results are 
presented in Figure l5T| where we see that the line segments have been substantially 
reduced. For a homogeneous object we can see what is going on also on the matrix 
level. The general form of our system matrix is A = I — GX, where I is a 37V x 37V 
identity matrix, X is the diagonal matrix, containing the grid values of the contrast 
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function, and G is a dense matrix produced by the integral operator with the Green's 
tensor kernel. When the contrast is homogeneous and equal to a scalar complex value 
X = e/eb — 1 inside the scatterer, we can write our matrix as A = I — \G. Now, let 
ip n be an eigenvector of G, i.e. Gip n = {5 n i\) n . Then, the eigenvalues of A are: 



An 



l-A, 



£ 
£b 



1 



and the eigenvalues of ArA are 

a:, = 



Eliminating j3 n we get 



A' — 



1 



«.(*-.) 



-(i-a„: 

£ 



A.-, 



(5.1) 



(5.2) 



(5.3) 



Wc have verified numerically that this transformation holds for the first two plots 
of Figure 14.21 (top) and Figure 15.11 Curiously, the seemingly quadratic amplification 
of the large eigenvalues in (|5.3[) is moderated by the inverse of the relative permit- 
tivity. That is why the regularized spectra are not as extended as one would fear. 
Unfortunately, such a simple proof is not possible with an inhomogeneous object, yet 
the effect of the regularization on the matrix spectrum in Figure l4~2l (top, right) and 
Figure [5J] (right) is obviously similar to the homogeneous cases. 

Table I5TT1 shows that such a regularized system is already better than the original 
one in terms of convergence. The relative speed-up, however, is limited to approx- 
imately 2.5 times due to the extra work involved in computing the matrix-vector 
products with A-&A. Moreover, no convergence improvement is observed with the 
negative permittivity object. 

To further accelerate the convergence and make the method work with the neg- 
ative permittivity objects as well (they are important in surface plasmon-polariton 
studies) we can now apply the deflation technique. For that we need to decide on the 
number of eigenvalues to be deflated. Obviously, it does not make sense to deflate 
beyond the first (from outside) dense cluster - deflation will only cost more time and 
memory, while the improvement will stagnate. After the regularization, the outer 
boundaries of the dense line segments have been all shifted inwards. However, due to 
nonlinear nature of the mapping (15.31) it is hard to tell exactly what an image of an 
arbitrary line segment would look like. Only for a lossless homogeneous scatterer are 
we able to predict that a real line segment [l,e/eb] will be mapped onto the real line 
segment [1, 1 + (1 — e/£b) 2 /(4e/£b)]j which contracts the original segment by approxi- 
mately four times. For a general line segment radially emerging from the point 1 + i0 
we have made a small numerical routine, which produces a discretized image according 
to (|5.3j) and finds the point with the largest absolute value - the radius of the circle 
beyond which all eigenvalues should be deflated. We have observed, by analyzing the 
spectra for various inhomogeneous objects, that the general regularization seems to 
transform each line segment according to (|5.3p . Therefore, in an inhomogeneous case, 
we can simply apply the aforementioned numerical mapping to each of the segments 
[1, e(x)/ef,], x G D, and choose the point with the maximum absolute value as the 
deflation radius. The deflation bounds obtained in this way are shown in Figure 15.11 
as solid red circles. 

The number of eigenvalues outside the deflation bound depends on the physics 
of the problem. For example, we know that at low frequencies (large wavelengths) 
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Fig. 5.2. Dependence of the number of eigenvalues to be deflated on the wavelength for three 
test scatterers (same medium parameters as in previous figures). Left: ratio ft of the number r of 
eigenvalues outside the deflation bound (circles with radius R, see Fiaure [5.1\l to the total number 
3-/V(A b ) of eigenvalues for a classical discretization rule (fixed number of points per smallest medium 
wavelength). Right: actual number r of eigenvalues to be deflated (with magnitude greater than the 
deflation bound R) and our analytical prediction based on an average /3 coefficient estimated from 
the left plot. 



all discrete eigenvalues are located within the convex hull of the essential spectrum 
[3J, and that at higher frequencies they tend to spread out. We have investigated this 
phenomenon in more detail and have found a pattern, similar to the one reported in 
[24] . As one can see from Figure IB~2l (left), the ratio of the number r of eigenvalues 
outside the deflation bound to the total number of eigenvalues 3N is roughly constant 
(for objects larger than At,), if at each wavelength we adjust the spatial discretization 
according to the classical rule (say, 15 points per smallest medium wavelength). Hence, 
with fixed medium parameters and a square computational domain with side a, the 
total number of unknowns and therefore eigenvalues grows as 3iV = 3[a 2 (a/Ab) 2 + 
2a(a/Ab) + 1] with a = 15 max 

W( £ '/ £ b)}- Thus, knowing f3 = r/(3N) for some 
wavelength, we can estimate r = 3Nf3 at any other wavelength. This empirical law 
comes very handy since solving a relatively low-frequency spectral problem can be 
really easy due to small system matrix dimensions. Otherwise, if there is no shortage 
of memory for the chosen wavelength, one can simply run the eigs routine increasing 
the number of recovered largest eigenvalues until the deflation bound obtained above 
is reached. 

Having a good estimate of the deflation parameter r we have applied the right 
preconditioner given by f|4. 1 2[) to the regularized system for the same three objects 
as in all previous examples. The largest eigenvalues/eigenvectors of the regularized 
system turn out to be very stable, so that one can safely apply the eigs routine 
with its tolerance set as high as 10 -2 . Though, we have not seen any substantial 
acceleration beyond the 10 -4 tolerance, which, in its turn, gave us a significant speed- 
up in the off-line work compared to the default machine precision tolerance. In all 
our experiments, the off-line time remained a small fraction of the total CPU time 
(for concrete numbers we refer to the example at the end of this Section). 

The convergence histories for a restarted GMRES with this preconditioner are 
presented in Figure [5T31 This figure, and the iteration counts in general, are slightly 
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Fig. 5.3. Convergence of restarted GMRES with the original Au = b (solid) and preconditioned 
ArAP~ 1 v = A^b (dashed) systems for the three physically distinct scatterers in the TE case. 
Red lines; e/eb = 16; Green lines: e/eb = — 16 + £1.5/ Blue lines: inhomogeneous object e/eb = 
(16,2.5 + i20). 



Table 5.1 

Restarted GMRES with fair memory usage for three different test scatterers, TE case. Original, 
deflated, regularized, and regularized-plus-delfated systems. 



System 


e/e b 


r, deflated 
eigenvalues 


restart 


CPU time, 
seconds 


Iterations, 
tol. 10~ 8 


Speed- 
up 


A 


16 





40 


350 


15258 


1.0 


AP- 1 


16 


6 


34 


306 


14002 


1.1 


A K A 


16 





40 


142 


4030 


2.5 


A R AP-i 


16 


14 


26 


117 


3626 


3 


A 


-16 + il.5 





40 


39 


1708 


1.0 


AP- L 


-16 + il.5 


7 


33 


36 


1659 


1.1 


A K A 


-16 + il.5 





40 


36 


1021 


1.1 


A R AP-" 


-16 + il.5 


30 


10 


16 


534 


2.4 


A 


16 and 2.5 + i20 





40 


275 


11823 


1.0 


AP- V 


16 and 2.5 + i20 


3 


37 


310 


13156 


0.9 


A R A 


16 and 2.5 + i20 





40 


107 


3032 


2.6 


A R AP-i 


16 and 2.5 + i20 


14 


26 


94 


2911 


3 



misleading when one deals with an iterative algorithm like restarted GMRES. For 
example, we see in Table 15.11 that, apart from the negative permittivity case, the 
iteration counts after deflation are only slightly better than what we have already 
achieved with the regularized system. Hence, one could conclude that the deflation is 
not worth the effort. However, the GMRES algorithm does not scale lineraly in time 
as a function of the restart parameter, i.e., while it may converge in less iterations 
for larger values of restart, the CPU time may not decrease as much. This has partly 
to do with the re-orthogonalization process, which takes more time for larger inner 
Krylov subspaces, and partly with the very bad spectral properties of our system 
matrix. 

Thus, to decide upon the benefits of deflation we compare execution times for var- 
ious values of restart parameter. In addition, we need to consider another important 
constraint - the available memory. Since a fair comparison can only be achieved if we 
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1 1.2 1.4 1.6 1.8 2 2.2 1 1.2 1.4 1.6 1.8 2 2.2 



a/A b a/A b 

Fig. 5.4. CPU time (left) and the achieved speed-up (right) as a function of wavelength (object 
size) for restart = r/2 (red, squares), restart = 2r (green, triangles), and restart = lOr (blue, 
circles). The permittivty of the scatterer is e/eb = 16. Solid lines - original system Au = b; dashed 
lines - preconditioned system A^AP _1 v = A^b. 



give the same amount of memory to both the original and the preconditioned systems, 
we give all the available memory to the GMRES algorithm in the former case, and split 
it between the restart and the deflation subspace of size r in the latter. Finally, we 
want to know if the efficiency of the proposed preconditioner depends on the physics 
of the problem, i.e., does it work better or worth for larger objects and contrasts? Of 
course, the system becomes larger and more difficult to solve for larger objects, and we 
have to be prepared to increase the value of restart for larger objects and contrasts 
to achieve any convergence at all. Hence, assuming that we have enough memory at 
least for deflation, we shall distinguish between the following three situations: limited 
memory, restart(j4Ry4P _1 ) < r; moderate memory, Testaj:t(A^AP^ 1 ) ~ r; large 
memory, restart(ARj4P _1 ) > r. 

We observe from the plots of Figure 15.41 that with limited memory, where we 
choose restart(ARv4P _1 ) = r/2 and restart(yl) = r, the original system initially 
outperforms the preconditioned one. As the object size grows, so does the size r of 
the deflation subspace. Since in that case the restart^RAP -1 ) parameter grows as 
well, at some point (for a/ A = 1.9 with this particular object) it reaches a value at 
which the preconditioned system breaks even and begins to outperform the original 
one for larger objects (maximum achieved speed-up is 2.4 times at a/Ab = 2.2). With 
moderate memory, restart(Ap;ylP _1 ) = 2r, restart(j4) = 3r, the preconditined 
system always performs better than the original one. Moreover, the relative speed-up 
Figure 15.41 (right) generally increases with the object size (maximum achieved speed- 
up is 5.1 times at a/Ab = 2, minimum - 2.1 times at a/Ab = !)• With large memory, 
i.e., for restart(ylR,AP _1 ) = 10> and restart(A) = llr, the speed-up may be as 
high as 46.6 times (a/Ab = 2.2) (minimum speed-up of 8.7 times at a/Ab = 1), and it 
also grows (on average) with the object size. 

With objects of finite extent changing the wavelength of the incident field is not 
the same as changing the object permittivity. Hence, we repeat the same numerical 
experiments as in Figure [5T5| but now for an object of fixed size, a/Ab = 1, and varying 
permittivity. The deflation parameter r is estimated using the same algorithm as 
above and turns out to grow only slightly as a function of permittivity, climbing from 
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Fig. 5.5. CPU time (left) and the achieved speed-up (right) as a function of permittivity for 
restart = r/2 (red, squares), restart = 2r (green, triangles), and restart = lOr (blue, circles). 
The scatterer is lossless with side length fixed at a/X^ = 1. Solid lines - original system Au = b; 
dashed lines - preconditioned system A^AP~ 1 v = A^b. 



13 to 17 for the permittivity changing between 10 and 20. At the same time the total 
number of unknowns as well as eigenvalues grows very fast. Therefore, deflation did 
not accelerate the iterative process here as much as when we changed the wavelength, 
and most of the speed-up must be due to the regularization. For this relatively 
small scatterer we see some improvement in convergence only with relatively large 
amounts of memory devoted to the restart parameter, i.e., restart(^4ri^4.P~ 1 ) > 2r. 
Moreover, preconditioned GMRES with rest&rt(A^AP^ 1 ) = r/2 did not converge 
for e/eb = 14 at all (we set this point arbitrarilly high in the plot). On the other 
hand, the maximum speed-up for restart = 2r was 4.1 times (ninimum 1.8 times), 
while for restart = lOr the maximum speed-up was 11.4 times (minimum 1.6 times). 
From the previous experiments we may expect that the speed-up will be larger with 
larger objects. 

As an example of a systematic application of the present preconditioner, we con- 
sider the problem, for which the DIE method is most suited, i.e., scattering on an 
object with continuously varying permittivity. The steps we suggest to follow are: 

1. Determine the total available memory: M bytes 

2. Fix the discretization rule at k points per smallest medium wavelength. 

3. Estimate f3 from low- frequency matrix spectra. 

4. Determine the maximal number of grid points N m!iX : 

• If a - size of the computational domain D, and N - number of grid 
points, then N = a 2 (a/Xb) 2 + 2a{a/\^) + 1, where 

a = fcmax{ Y / (£'(x)/eb)}, x e D. 

• Memory needed for the system, right-hand side, and the unknown vector: 
Ma = 720N bytes (one complex number - 16 bytes). 

• Memory required for deflation and GMRES: 

Afrroc = 16r 2 + 48Ar + 48iVrestart bytes, where r = 3N/3 and 
restart = xr, i.e., M Prcc = 144(/3 2 + (3 + xf3)N 2 bytes. 

• The largest affordable grid follows from Ma + Mp ICC ps M: 
A^max w [V900 + M(/3 2 + p + xj) - 30]/[12(/3 2 + + xj3)}. 

5. Determine the maximum affordable object size (smallest wavelength) via 
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xi xi 

Fig. 5.6. Scattering of a TE plane wave on an inhoraogeneous object with continuously varying 
permittivity. Left - permittivity profile (side length a/Ab = i). Right - computed electric field 
intensity \E\\ 2 + |i?2| 2 . Preconditioning gives 14 times speed-up for this problem. 

[a/A b ]max = (V^max - l)/«- 

The scatterer is depicted in Figure [576] (left). The relative permittivity function profile 
is given by 

e(x)/eb = 10 + 5 sin(47ra;i/a) sin(47ra;2/a), (5.4) 

where the origin of the coordinate system is in the middle of the computational do- 
main. We calculated j3 for this scatterer to be approximately 0.0011. Having M — 
4 GB of memory available and choosing k = 15 and restart(Aft^4P _1 ) = 8r, from 
the above algorithm we determine the largest affordable size with this permittivity as 
[a/Ab]max = 4. Hence, the total number of unknowns is 3N = 162867 and the defla- 
tion parameter is r — 180. The GMRES was restarted at restart(,4 R AP -1 ) = 1440 
and restart(A) = 1620. The results of this test are as follows: original system - 
162692 seconds of CPU time, preconditioned system - 11620 seconds. The speed-up 
in pure 'on-line' calculations is 14 times, as we would expect from our previous ex- 
periments with the homogeneous test object. The off-line work with the tolerance of 
eigs set to 10~ 4 took only 775 seconds, so that even with this time included we get 
a 13.1 times speed-up. The computed electric field is shown in Figure [5761 (right). 

The CPU timing is provided here merely for illustration and should not be used 
for direct comparison against the local methods, which was not the purpose of this 
paper either. For the particular problems considered above one could, for instance, 
substantially accelerate the calculations and free some valuable memory by computing 
the electric field only. On the other hand, the presented CPU timing data give us 
a good relative estimate of what to expect with objects having both electric and 
magnetic contrasts so that both electric and magnetic fields must indeed be computed 
at the same time. 

Finally we should mention that, while we have been focusing here on the restarted 
GMRES, obviously, any other memory-efficient solver, which works with the original 
system Au = b, could also be applied to the preconditioned system AnAP~ 1 v — A^b. 
We have considered, for instance, Bi-CGSTAB, restarted GCR, and IDR. Sometimes 
we could achieve faster convergence by tuning the restart method of GCR, although, 
this did not happen with all types of scatterers (only with a homogeneous lossless 
one). Moreover, with some scatterers (inhomogeneous with losses) GCR and other 
algorithms would simply stagnate. Thus, for the moment, the restarted GMRES 
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remains the most robust solver for the present problem, whereas other algorithms 
require more tuning and a more systematic study. 

6. Conclusions. The analysis of the domain integral operator of the transverse 
electric scattering presented here helps to understand the reasons behind the ex- 
tremely slow convergence of iterative methods often observed with this method. Un- 
like the TM case, which is equivalent to the Helmholtz equation and can be easily 
accelerated by deflating the largest eigenvalues, the TE case showed no improvement 
with this type of preconditioncr. Direct numerical comparison of the TE and TM 
matrix spectra indicates that the nontrivial essential spectrum of the TE operator 
results in extended and very dense eigenvalue clusters, that cause stagnation of the 
deflation process with high-contrast objects. We have derived an analytical multi- 
plicative rcgularizcr and demonstrated on various test scatterers that its discretized 
version robustly reduces the extent of the troublesome dense clusters. One can view 
the obtained regularizer as a generalization of the Calderon identity, previously ap- 
plied in the preconditioning of boundary integral equations. Notably, the discretized 
regularizer does not require any extra computer memory and its action on a vector 
can be computed at FFT speed. Applying deflation of the largest-magnitude eigen- 
values on the regularized system we could achieve up to 46.6 times acceleration of 
the restarted GMRES with relatively large memory, and up to 5.1 times with modest 
memory. The off-line work typically takes only a small fraction of the total time, 
since one can apply a very rough tolerance in the eigs algorithm when recovering 
the largest- magnitude eigenvalues/eigenvectors. Somewhat surprisingly, this precon- 
ditioner showed a tendency to become more efficient (on average) with the increase 
in the object size and permittivity. 
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